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1.  Introduction 


Three-dimensional  (3-D)  operations,  such  as  rotations,  translations,  and  intersections,  are  tools 
that  are  essential  for  many  types  of  scientific  modeling.  However,  the  C++  programming 
language  does  not  natively  perfonn  them.  This  report  documents  a  set  of  functions,  written  in 
C++,  that  can  be  used  to  perform  3-D  rotations,  translations,  and  intersections.  All  of  the 
functions  have  been  grouped  together  in  the  namespace  y3DOps.  A  summary  sheet  that 
reproduces  the  entire  y3DOps  namespace  is  located  at  the  end  of  this  report. 


2.  Translation  of  a  Point  in  Space 


2.1  Derivation 

Let  the  position  vector  p  represent  an  arbitrary  point  in  space,  where 

P  =  Pxx  +  Pyy  +  pJ.  (1) 

Furthermore,  let  d  represent  a  displacement  vector,  where 

d  =dxx  +  dyy  +  dzz .  (2) 

If  p  is  used  to  represent  p  after  it  has  been  translated,  then 

P  =  (px  +  dx  )x  +  (py  +  dy  )y  +  (p:  +  dz  )z .  (3) 


2.2  C++  Implementation 

TranslateQ  Code 


inline  void  Translate(//<===============================PERFORMS  A  TRANSLATION 

double  p[3 \J/< . -COORDINATES  TO  TRANSLATE  (MODIFIED  BY  THIS  FUNCTION) 

const  double  d  [3]  ){//< . DISPLACEMENT  VECTOR 

p[0]+=d[0]  ,  p[l]+=d[l]  ,  p[2]+=d[2] ; 

}//  YAGENAUT@GMAIL.COM - LAST~UPDATED~02MAY2013 - 


Translate()  Parameters 

p  p  is  a  three-element  array  that  stores  the  position  vector  that  is  described  by  equation  1 . 
Note  that  p  is  modified  by  the  Translate()  function,  as  described  by  equation  3. 

d  d  is  a  three-element  array  that  stores  the  displacement  vector  that  is  described  by 
equation  2.  d  determines  the  amount  and  direction  by  which  p  is  translated. 
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Translate!)  Example 

Figure  1  shows  point  p  being  translated  to  a  new  position  ( p' )  by  displacement  vector  d  . 


/v 

Z 


Figure  1.  Translate!)  example. 

Let  p  =  {l,l,0}  and  d  =  {- 1,-1, 2} .  Point  /)'  can  be  found  by  using  the  Translate!)  function,  as 
shown  in  the  following  sample  code. 


#include  <cstdio>// . printfQ 

#include  "y_3d_ops . h"// . Translate! ) 

int  main(){ 

double  p[3]={l,l,0}; 


double  d[3]={-l, -1,2}; 
y3D0ps: :Translate(p,d); 

printf("p[0]=%f,  p[l]=%f,  p[2]=%f\n", p[0] j p[l] , p[2] ) ; 

}//  VAGENAUT@GMAIL.COM - LAST~UPDATED~02MAY2013 - 


OUTPUT: 

p[0]=0.000000j  p[l]=0. 000000,  p [ 2 ] =2 . 000000 


3.  Rotation  of  a  Point  About  an  Arbitrarily  Positioned  Axis 


3.1  Derivation 

Suppose  that  the  unit  vector  v  is  used  to  define  an  arbitrary  axis  about  which  a  point  in  space 
will  be  rotated. 
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v  =  vxx  +  vyy  +  vzz. 


(4) 


Rodrigues’  rotation  fonnula  can  be  used  to  construct  a  rotation  matrix,1  R  ,  that  can  be  used  to 
perform  a  rotation  about  v  by  an  angle  6 .  Note  that  the  direction  of  the  rotation  will  be 
determined  by  the  right-hand-thumb  rule  (when  the  right  thumb  is  pointed  in  the  direction  of  v , 
the  curled  lingers  of  the  right  hand  will  point  in  the  direction  of  the  rotation). 


VxO-cJ  +  Ce 

vxvy(l-co)~v*se 

VxV:(\-Cg)+VySg 

VxVy(l-Cd)+VzS0 

V;(l~Cg)+Cg 

(5) 

VxVz^~ce)-Vyse 

vyvAl-ce)  +  vxsa 

V2z(\-Cg)+Cg 

where 

cg  =  cos(<9)  (6) 

and 

sg  =  sin(<9) .  (7) 

Let  the  position  vector  p  locate  an  arbitrary  point  in  space,  where 

p  =  px* + pyy + pJ  •  (8) 

Let  the  position  vector  o  locate  the  origin  of  the  rotation  axis  defined  by  v . 

o=oxx  +  oyy  +  o2z.  (9) 

The  translation-rotation-translation  sequence  described  by  equation  10  can  be  used  to  find  p' , 
where  p'  is  used  to  represent  p  after  it  has  been  rotated  about  v . 

p  =  R(p  —  o)+o  .  (10) 

3.2  C++  Implementation 

Two  functions  are  used  to  perfonn  3-D  rotations.  The  first  function,  RMatrix(),  calculates  the 
rotation  matrix  that  is  presented  in  equation  5.  The  second  function,  Rotate(),  perfonns  the 
rotation  that  is  presented  in  equation  10.  Breaking  the  calculation  into  two  functions  allows 
functions  that  rotate  objects  containing  more  than  one  point  to  be  written  in  a  manner  that 
doesn’t  sacrifice  performance. 

A  special  note  about  rotations  -  Although  the  rotation  functions  presented  in  this  report  were 
designed  to  be  used  to  rotate  position  vectors,  they  can  also  be  used  with  other  types  of  vectors. 
For  example,  suppose  that  a  multi-point,  rigid  object  has  a  velocity  vector  associated  with  it. 

'Mason,  M.  T.  Mechanics  of  Robotic  Manipulation',  Massachusetts  Institute  of  Technology:  United  States,  2001,  (page  46, 
equation  3.26). 
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That  object  can  be  rotated  by  first  creating  a  rotation  matrix  using  the  RMatrix()  function.  That 
rotation  matrix  can  then  be  supplied  to  the  Rotate()  function  to  rotate  all  of  the  object’s  points. 
Using  the  same  rotation  matrix,  the  Rotate()  function  can  also  be  used  to  rotate  the  object’s 
velocity  vector.  However,  the  origin  (i.e.,  the  second  argument  of  the  Rotate()  function)  will 
need  to  be  set  to  {0,0,0}  for  the  call  that  rotates  the  velocity  vector.  Other  types  of  vector 
quantities,  such  as  torque,  angular  velocity,  and  rotation  axes  can  similarly  be  rotated.  However, 
the  origin  will  typically  need  to  be  zeroed  for  vectors  that  are  not  displacement  vectors. 

RMatrix})  Code 


inline  void  RMatrix(//<========================CALCULATES  A  3D  ROTATION  MATRIX 

double  R[9],//< . ROTATION  MATRIX  (CALCULATED  BY  THIS  FUNCTION) 

const  double  v[3],//< . THE  AXIS  OF  ROTATION  (V  MUST  BE  A  UNIT  VECTOR) 

double  rads){//<--THE  ANGLE  OF  ROTATION  (CCW  ABOUT  V,  USE  RIGHT-HAND  RULE) 
double  c=cos ( rads ), d=l-c, s=sin (rads), v0=v[0] ,vl=v[l] ,v2=v[2] ; 

R[0]=v0*v0*d+  c  ,  R[l]=v0*vl*d-v2*s  ,  R[2]=v0*v2*d+vl*s  ; 

R[3]=vl*v0*d+v2*s  ,  R[4]=vl*vl*d+  c  ,  R[5]=vl*v2*d-v0*s  ; 

R[6]=v2*v0*d- vl*s  ,  R[7]=v2*vl*d+v0*s  ,  R[8]=v2*v2*d+  c  ; 

}// —  ~~YAGENAUT (5)GMAI L . LAST ~U  P  DAT  E  D~0  2MAY 201 


RMatrix()  Parameters 

R  R  is  a  nine-element  array  that  stores  the  rotation  matrix  that  is  described  by  equation  5. 
Note  that  R  is  modified  by  the  RMatrix()  function.  R  is  intended  to  be  used  as  the  third 
argument  of  the  Rotate})  function. 

v  v  is  a  three-element  array  that  stores  the  rotation  axis  that  is  described  by  equation  4. 

Note  that  v  must  be  a  unit  vector. 

rads  rads  is  used  to  represent  the  angle  (in  radians)  of  the  rotation.  As  described  in  section 
3.1,  the  direction  of  the  rotation  is  determined  by  the  right-hand-thumb  rule. 

Rotate})  Code 


inline  void  Rotate(//<=====================================PERFORMS  A  ROTATION 

double  p [3 ],//< . COORDINATES  TO  ROTATE  (MODIFIED  BY  THIS  FUNCTION) 

const  double  o[3],//< . -THE  ORIGIN  OF  THE  AXIS  OF  ROTATION 

const  double  R[9]){//< . A  ROTATION  MATRIX  (FROM  RMatrixQ) 

double  t0=p[0]-o[0],tl=p[l]-o[l],t2=p[2]-o[2];// . translate  to  the  origin 

p[0]=R[0]*t0+R[l] *tl+R[2] *t2+o[0] ; // . rotate,  then 

p[l]=R[3]*t0+R[4] *tl+R[ 5] *t2+o[l] ; //  reverse 

p[2]=R[6]*t0+R[7] *tl+R[8] *t2+o[2] ; //  translation 

/  j  ‘'"'"'"N'Y  AG  E  N  AU  T  (S)GMA  XL*  C  0^1  ^  r^/  ru  nu  ru  r^/r^t  nu  ru  nu  nu  ru  nu  ru  nunu  ru  nu  ru  nu  ru  nu  ru  LAST  ^  LJ  P  DAT  E  D  ^  0  2  MAY  2013  iV'V'N/'V'N/'V 


Rotate})  Parameters 

p  p  is  a  three-element  array  that  stores  the  position  vector  that  is  described  by  equation  8. 
Note  that  p  is  modified  by  the  Rotate})  function,  as  described  by  equation  10. 
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o  o  is  a  three-element  array  that  stores  the  position  vector  that  is  described  by  equation  9. 
o  is  the  origin  of  v,  the  axis  about  which  p  is  rotated. 

R  R  is  a  rotation  matrix  that  has  been  precalculated  using  the  RMatrix()  function. 

Rotate()  Example 

Figure  2  shows  point  p  being  rotated  about  the  axis  v  to  a  new  position  (/?'). 


/v 

Z 


Figure  2.  Rotate()  example. 

Let  p  =  {2,0,0}  and  v  =  {0,0,l} .  Furthermore,  let  the  origin  of  v  be  at  o  =  {0,0,0} ,  and  the  angle 
of  rotation  be  n  12.  Point  p'  can  be  found  by  first  using  the  RMatrix()  function  to  calculate  a 
rotation  matrix,  then  using  the  RotateQ  function  to  perfonn  the  rotation. 


#include  <cstdio>// . printfQ 

#include  "y_3d_ops .  h "// . RMatrix( ) , Rotate ( ) 

int  main(){ 

double  p[3]={2j0j0}; 

double  v[3]={0j0,l};// . note  v  must  be  a  unit  vector 

double  o[3]={0j 0j 0}; 

double  R [ 9 ] ; /*< - */y3D0ps : : RMatrix(Rj v, 3 . 14159265358979/2) ; II...  rotation  matrix 
y3D0ps : : Rotate(pj  o, R) ; 

printf ("p[0]=%f,  p[l]=%f,  p[2]=%f\n", p[0] j p[l] , p[2] ) ; 

}//—  AG  E  N  AU  T  ^)GF1A  XL*  C  0M,v,v,v,vrvj,vrv,v,>j,N/,v,Ny,v,v,vrv,v,N^,v,N/rw,v,v,v,>j,v,v  LAS  T  ^  LJ  P  DAT  ED^G  ^  MAY  2013 


OUTPUT: 

p[0]=0.000000j  p [ 1 ] =  2 . 000000,  p[2]=0. 000000 
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4.  Intersection  Between  a  Line  and  a  Plane 


4.1  Derivation 

Suppose  that  a  line  S  passes  through  the  points  S0  and  5, ,  where 

S0  =  S0,xx  +  S0yy  +  S0J  and  S,  =  Sl  xx  +  SUyy  +  Suz  .  (11) 

Let  ps  represent  a  point  that  lies  on  S  .  S0  and  Sl  can  be  used  to  construct  a  parametric 
equation  for  ps  as  a  function  of  the  parameter  t0 : 

ft=s,+(s,-s0)r0.  (12) 

The  parameter  t0  represents  the  scaled  distance  from  S0  to  Sl  along  S  .  Thus  if  t0  =  0 ,  ps  is 
located  at  S0 .  If  t()  =  1 ,  ps  is  located  at  Sr . 

Next,  suppose  that  a  plane  T  passes  through  the  points  T0,  Tx,  and  T2 ,  where 


To=ToJ  +  ToJ  +  ToJ 


Ti  =  TiJ  +  TiJ  +  TiJ ,  and  T2  =  Tlxx  +  T2  y  +  T2J  . 


Let  pT  represent  a  point  that  lies  on  f.  T0 ,  Ij ,  and  T2 ,  can  be  used  to  construct  a  parametric 
equation  for  pT  as  a  function  of  the  parameters  t]  and  t2 : 

Pt  =  T0  +  (^  +{f2  -T0)t2 .  (14) 


Similar  to  t0 ,  t]  and  t2  represent  scaled  distances  from  T0  to  Tx  and  from  T0  to  T2 ,  respectively. 
The  point  of  intersection  between  S  and  T  occurs  where  ps  is  equal  to  pT .  Thus, 

•So  +(s.  -S„k  =f«  +(f,  -f  h  +(f,  (15) 

Rearranging  terms, 

-f0  =  (s0 -Sjt 0  +(ft  -f0)t1+(f2  -f0)t2.  (16) 

This  can  be  written  in  matrix  form  as 
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Solving  for  t  , 


*0 

e  _e  T  —T  T  —T 

°0,x  °l,x  1l,x  10,x  1 2,x  10,x 

-1 

1 

E~?" 

1 

* 

o*~ 

Co 

1 _ 

h 

= 

1 

i 

crT 

1 

o' 

So,y-To, 

^2  _ 

i 

N 

1 

E*? 

1 

N. 

of 

1 

O*' 

Co 
_ 1 

\.~T0,Z_ 

(18) 


Recall  that  the  vector  t  contains  the  parameters  from  equations  12  and  14.  Thus,  once  t  is 
known,  t0  can  be  substituted  into  equation  12  to  find  the  point  of  intersection  between  S  and  T. 

Note  that  if  the  three-by-three  matrix  defined  in  equation  18  is  noninvertible,  then  S  is  parallel 
to  T. 


Equation  14  can  be  used  to  find  a  set  of  conditions  that  test  whether  S  intersects  T  inside  or 
outside  the  triangle  defined  by  the  three  points  T0 ,  T[ ,  and  T2 .  Note  the  colored  bands  shown  in 
figure  3.  The  blue  band  represents  the  region  where  0<h<l  ,  and  the  yellow  band  represents 
the  region  where  0<t2<l. 


The  region  where  the  two  bands  overlap  is  where  both  sets  of  conditions  are  met.  Thus,  the 
conditions  presented  in  equation  19  can  be  used  to  determine  whether  or  not  line  S  intersects  the 
plane  inside  the  parallelogram  shown  in  green  in  figure  3. 

0  <  tx  <  1  and  0  <  t2  <  1 .  (19) 

The  additional  test  shown  in  equation  20  can  be  used  to  determine  whether  or  not  line  S 
intersects  the  plane  inside  the  triangle  defined  by  the  three  points  T0,  T[ ,  and  T2 . 
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(fi+Q)<  I-  (20) 

Finally,  if  the  condition  presented  in  equation  2 1  is  met,  then  S  intersects  T  between  S0  and  S1 . 

0<t0<l.  (21) 


4.2  C++  Implementation 

Two  functions  are  used  to  find  line-plane  intersections.  The  first  function,  IParameters(), 
calculates  a  three-element  array  that  is  the  solution  to  equation  18.  The  second  function, 
Intersect(),  calculates  the  point  of  intersection  between  a  line  and  a  plane. 

Because  there  is  a  chance  that  the  three-by-three  matrix  shown  in  equation  1 8  will  be  singular,  a 
Boolean  that  indicates  whether  or  not  a  solution  is  valid  is  returned  by  the  IParameters() 
function.  The  intersection  tests  described  in  equation  19  are  perfonned  by  the  Intersect() 
function. 

Note  that  the  Intersect()  function  tests  for  intersection  between  a  parallelogram  and  a  line.  If  a 
test  for  intersection  between  a  triangle  and  a  line  is  desired,  then  the  additional  condition 
presented  in  equation  20  must  be  tested  for.  If  a  test  for  intersection  between  a  triangle  and  a  line 
segment  is  desired,  then  the  additional  condition  presented  in  equation  2 1  must  be  tested  for. 
Neither  test  is  performed  by  the  Intersect()  function. 

IParametersQ  Code 


inline  bool  IParameters(//<=============PARAMETERS  FOR  LINE-PLANE  INTERSECTION 

double  t [3 ],//< . INTERSECTION  PARAMETERS  (CALCULATED  BY  THIS  FUNCTION) 

const  double  T[9],//<--A  TRIANGLE  {T0XJT0YJT0ZJT1XJT1YJT1ZJT2XJT2YJT2Z} 

const  double  S[6],//< . A  LINE  {SOX.SOY.SOZ.SIX.SIY^SIZ} 

double  e=lE-9){//< - CUTOFF  VALUE  FOR  DETERMINING  IF  S  &  T  ARE  PARALLEL 

double  A0=S[0] - S [ 3 ]  ,  A1=T[3]-T[0]  ,  A2=T[6] -T[0] ,// . A  in  y=A*t 

A3=S[1] -S[4]  ,  A4=T[4] - T [ 1 ]  ,  A5=T[7] -T[l] , 

A6=S[2] - S [ 5 ]  ,  A7=T[5] -T[2]  ,  A8=T[8] -T[2] ; 
double  d=A0*A4*A8-A0*A5*A7  +  A1*A5*A6-A1*A3*A8  +  A2*A3*A7-A2*A4*A6; // . . d= | A | 
if (fabs(d)<e)return  false; //.... =>  S  &  T  are  parallel  (and  t  is  meaningless) 
double  B0=(A4*A8-A5*A7)/d , Bl=(A2*A7-Al*A8)/d , B2=(Al*A5-A2*A4)/d , / / . A  inverse 
B3=(A5*A6-A3*A8)/d , B4=(A0*A8-A2*A6)/d , B5=(A2*A3-A0*A5)/d , 
B6=(A3*A7-A4*A6)/d j  B7=(Al*A6-A0*A7)/d , B8=(A0*A4-Al*A3)/d; 

double  y0=S[0]-T[0] ,yl=S[l]-T[l] ,y2=S[2]-T[2];// . y  in  y=A*t 

t [0] =B0*y0+Bl*yl+B2*y2  ,  t[l]=B3*y0+B4*yl+B5*y2  ,  t[2]=B6*y0+B7*yl+B8*y2; 

return  true;// . =  >  S  &  T  intersect 

"^11  •v'v'v'N'YAG  ENAUT^IGMAI  L  •  LAST  ,s'UPDATED,''02F1AY2013‘v,n,,v,,n','<,v 


IParameters()  Parameters 

t  t  is  a  three-element  array  that  stores  the  parameters  described  in  equation  18.  Note  that 
t  is  modified  by  the  IParametersQ  function. 
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T  is  a  nine-element  array  that  stores  the  triangle  that  is  defined  by  equation  13. 
T=  {To?  ,T0v,  T0  ,  T.  ,T,  ,  Tx  ,  T2  ,  T2  ,  T2  }  . 


S  is  a  six-element  array  that  stores  the  line  that  is  defined  by  equation  1 1 . 

S=  {Sq,x  ’Sq  ,y>  *^0,z  !  j  S,  ,  S,z  }  . 


e  e  is  the  cutoff  value  for  testing  whether  or  not  S  and  T  are  parallel.  If  the  detenninant 
of  the  matrix  in  equation  18  is  less  than  e,  then  S  and  T  are  considered  to  be  parallel. 
The  default  value  of  e  is  10  9 . 


IParameters()  Return  Value 

IParameters()  returns  false  if  S  is  parallel  to  T.  A  return  value  of  false  indicates  that  t  has  not 
been  calculated  and,  thus,  should  not  be  passed  to  the  Intersect()  function. 

Intersect()  Code 


inline  bool  Intersect(//<=============CALCULATES  LINE-PLANE  INTERSECTION  POINT 

double  x[3],//< . POINT  OF  INTERSECTION  (CALCULATED  BY  THIS  FUNCTION) 

const  double  t[3],//< . -INTERSECTION  PARAMETERS  (FROM  IParameters( ) ) 

const  double  S[6] >{//< . A  LINE  {SOX.SOY.SOZ.SIX.SIY^SIZ} 

for(int  i=0; i<3;++i)x[i]=S[i]+(S[i+3] - S [ 1 ] )*t[0] ; 

return  t[l]>0&&t[l]<l&&t[2]>0&&t[2]<l;// . S  intersects  parallelogram? 

}//  YAGENAUT@GMAIL  .  COM - LAST~UPDATED~02MAY2013 - 


Intersect()  Parameters 

x  x  is  a  three-element  array.  Note  that  a  purpose  of  the  function  is  to  calculate  x. 

t  t  is  a  parameter  list  that  has  been  precalculated  using  the  IParameters()  function. 

S  S  is  a  six-element  array  that  stores  the  line  that  is  defined  by  equation  1 1 . 


Intersect()  Return  Value 

Intersect()  returns  true  if  line  S  intersects  the  parallelogram  shown  in  figure  3.  Note  that 
additional  checks  can  be  performed  to  determine  whether  or  not  the  point  of  intersection,  x  lies 
within  the  triangle  defined  by  T[0],  T[l],  and  T[2],  or  between  the  points  on  line  S,  S[0]  and 
S[l].  See  equations  20  and  2 1  for  the  additional  checks. 

Intersect^)  Example 

Figure  4  shows  the  point,  p  ,  where  line-segment  S  intersects  triangle  T . 
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Figure  4.  IntersectQ  example. 

Let  S0  =  2},  =  {l,l,2},  L0  =  {0,0,0},  7j  =  {4,0,0},  and  =  {0,4,0}.  Point  ]3  can  be  found 

by  first  calling  the  IParameters()  function,  then  using  the  result  in  the  Intersect()  function. 

#include  <cstdio>// . printfQ 

#include  "y_3d_ops . h"// . IParameters( ) , Intersect ( ) 

int  main(){ 

double  S[6]={ljlj  -2, 1,1, 2};// . a  line  segment 

double  T[9]={0J0J0J4J0J0J0J4J0};// . a  triangle 

double  t[3];/*<-*/y3DOps: :IParameters(t,T,S);// . intersection  parameters 

double  p[3];/*<-*/y3D0ps: :Intersect(p,t,S);// . point  of  intersection 

printf("p[0]=%f ,  p[l]=%f,  p[2]=%f\n",p[0],p[l],p[2]); 

}// —  s Y  AG  E  N  AU  T  ^)GMA  XL*  C  LAS  T  ^  LJ  P  DAT  ED^G  2  MAY  2013  'vrw,v'v,v'v 


OUTPUT: 

p[0]=l. 000000,  p[l]=l. 000000,  p[2]=0. 000000 
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5.  Summary 


A  summary  sheet  is  presented  at  the  end  of  this  report.  It  presents  the  y3DOps  namespace, 
which  contains  the  five  functions  that  are  described  in  detail  in  sections  2,  3,  and  4  of  this  report. 
Also  presented  are  two  examples  that  demonstrate  the  versatility  of  the  functions  described  in 
this  report.  The  first  uses  the  Rotate()  function  to  calculate  a  set  of  points  that  defines  the  surface 
of  a  sphere.  The  second  uses  the  Intersect()  function  to  shade  the  blades  of  a  pinwheel.  Both 
functions  create  text  files  that  contain  all  of  the  information  needed  to  create  the  two  images 
presented  in  the  summary  sheet.  The  actual  images  that  are  presented  were  created  using  the 
yGraphics  namespace,  which  will  be  documented  in  a  future  report. 
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y3DOps  Summary 


y_3d_ops . h 


ttifndef  Y_3D_0PS_H_ 

#define  Y_3D_0PS_H_ 

#include  <cmath>// . sin(),cos(),abs() 

namespace  y3D0ps{// - ALL  FUNCTIONS  ASSUME  RIGHT-HANDED  CARTESIAN  COORDINATES 


inline  void  Translate(//<===============================PERFORMS  A  TRANSLATION 

double  p[3],//< . -COORDINATES  TO  TRANSLATE  (MODIFIED  BY  THIS  FUNCTION) 

const  double  d[3]  ){//< . . . . . DISPLACEMENT  VECTOR 

p[0]+=d[0]  ,  p[l]+=d[l]  ,  p[2]+=d[2] ; 

}// - YAGENAUT@GMAIL.COM - LAST~UPDATED~02MAY2013 - 

inline  void  Rotate(//<=====================================PERFORMS  A  ROTATION 

double  p [ 3 ] , //< - COORDINATES  TO  ROTATE  (MODIFIED  BY  THIS  FUNCTION) 

const  double  o[3],//< . . THE  ORIGIN  OF  THE  AXIS  OF  ROTATION 

const  double  R[9]){//< - - A  ROTATION  MATRIX  (FROM  RMatrixQ) 

double  t0=p[0]-o[0],tl=p[l]-o[l],t2=p[2]-o[2];// . translate  to  the  origin 

p[0]=R[0]*t0+R[l]*tl+R[2]*t2+o[0] ;// . rotate,  then 

p[l]=R[3]*t0+R[4]*tl+R[5]*t2+o[l] ;//  reverse 

p[2]=R[6]*t0+R[7]*tl+R[8] *t2+o[2] ;//  translation 

}// - YAGENAUT@GMAIL.COM - LAST~UPDATED~02MAY2013 - 

inline  void  RMatrix(//<========================CALCULATES  A  3D  ROTATION  MATRIX 

double  R [9 ] ,  //< - ROTATION  MATRIX  (CALCULATED  BY  THIS  FUNCTION) 

const  double  v[3],//< - THE  AXIS  OF  ROTATION  (V  MUST  BE  A  UNIT  VECTOR) 

double  rads){//<- -THE  ANGLE  OF  ROTATION  (CCW  ABOUT  V,  USE  RIGHT-HAND  RULE) 
double  c=cos(rads),d=l-c, s=sin(rads),v0=v[0] , vl=v[l] , v2=v[2] ; 

R[0]=v0*v0*d+  c  ,  R[l]=v0*vl*d-v2*s  ,  R[2]=v0*v2*d+vl*s  ; 

R[3]=vl*v0*d+v2*s  ,  R[4]=vl*vl*d+  c  ,  R[5]=vl*v2*d-v0*s  ; 

R[6]=v2*v0*d-vl*s  ,  R [7 ]=v2*vl*d+v0*s  ,  R[8]=v2*v2*d+  c  ; 


}// - YAGENAUT@GMAIL.COM - LAST-UPDATED-02MAY2013 - 

inline  bool  Intersect (//<=============CALCULATES  LINE-PLANE  INTERSECTION  POINT 

double  x[3],//< . POINT  OF  INTERSECTION  (CALCULATED  BY  THIS  FUNCTION) 

const  double  t[3],//< . . INTERSECTION  PARAMETERS  (FROM  IParameters () ) 

const  double  S[6]){//< . . A  LINE  {S0X,S0Y,S0Z,S1X,S1Y,S1Z> 

for(int  i=0;i<3;++i)x[i]=S[i]+(S[i+3] - S [ i ] )*t[0]; 

return  t[l]>0&&t[l]<l&&t[2]>0&&t[2]<l;// . S  intersects  parallelogram? 

}// - YAGENAUT@GMAIL.COM - LAST-UPDATED-02MAY2013 - 

inline  bool  iIParamete^(//<=============PARAMETERS  FOR  LINE-PLANE  INTERSECTION 

double  t[3],//< . INTERSECTION  PARAMETERS  (CALCULATED  BY  THIS  FUNCTION) 

const  double  T[9],//<--A  TRIANGLE  {T0X,T0Y,T0Z,T1X,T1Y,T1Z,T2X,T2Y,T2Z> 

const  double  S[6],//< . . A  LINE  {S0X,S0Y,S0Z,S1X,S1Y,S1Z} 

double  e=lE-9){//< - CUTOFF  VALUE  FOR  DETERMINING  IF  S  &  T  ARE  PARALLEL 

double  A0=S[0]-S[3]  ,  A1=T[3]-T[0]  ,  A2=T[6]-T[0],// . A  in  y=A*t 

A3=S[1]-S[4]  ,  A4=T[4] -T[l]  ,  A5=T[7] -T[l] , 

A6=S[2]-S[5]  ,  A7=T[5]-T[2]  ,  A8=T[8] -T[2] ; 
double  d=A0*A4*A8 - A0* A5 *A7  +  A1*A5*A6-A1*A3*A8  +  A2*A3*A7-A2*A4*A6; // . . d= | A | 
if (fabs(d)<e)return  false;//. .. .=>  S  &  T  are  parallel  (and  t  is  meaningless) 
double  B0=(A4*A8-A5*A7)/d,Bl=(A2*A7-Al*A8)/d,B2=(Al*A5-A2*A4)/d,//.A  inverse 
B3=(A5*A6-A3*A8)/d,B4=(A0*A8-A2*A6)/d,B5=(A2*A3-A0*A5)/d, 
B6=(A3*A7-A4*A6)/d,B7=(Al*A6-A0*A7)/d,B8=(A0*A4-Al*A3)/d; 

double  y0=S[0] -T[0] ,yl=S[l] -T[l],y2=S[2] -T[2] ;// . y  in  y=A*t 

t [0]=B0*y0+Bl*yl+B2*y2  ,  t [l]=B3*y0+B4*yl+B5*y2  ,  t [2]=B6*y0+B7*yl+B8*y2; 

return  true;// . =>  S  &  T  intersect 

}// - YAGENAUT@GMAIL.COM - LAST-UPDATED-02MAY2013 - 


IMAGE  1 

image  created  from  sphere.txt. 


IMAGE  2 

image  created  from  pinwheel.txt. 


EXAMPLE  1 


#include  <algorithm>// . copy() 

#include  <cstdio>// . FILE, f reopen (), stdout,printf( ),fclose( ) 

#include  "y_3d_ops .  h"// . RMatrix( ) ,  Rotate( ) 

int  main( ){//<====================================USE  ROTATIONS  TO  DRAW  A  SPHERE 

using  namespace  y3D0ps; 
using  std : : copy; 

double  x[3]={l,0,0}  ,  y[3]={0,l,0}  ,  z[3]={0,0, 1}  ,  o[3]={0,0,0},a[3]={0,l,0}; 

double  H[200][23][3]={0,0,1>  ,  V[200] [23] [3]={1,0,0};// . lat  &  long  lines 

double  Ry[9];/*<-*/RMatrix(Ry,y,3.14159/(23+l)); 
double  Rz [9] ;/*<-*/RMatrix(Rz,z, 2*3. 14159/ (200-1)); 
double  R3[9];/*<-*/RMatrix(R3,z,3.14159/(23)); 
for(int  j=0; j<23;++j){ 

if ( j !=0)  Rotate (copy (H[0] [j-l]jH[0][j-l]+3,H[0][j])-3,o,Ry); 
for(int  i=l;i<200;++i)Rotate(copy(H[i-l] [j] ,H[i-l] [j]+3,H[i] [j])-3,o,Rz); 
if ( j !=0)  Rotate ( copy ( V[ 0] [j-l],V[0][j-l]+3,V[0][j])-3,o, R3), Rotate ( a, o,R3); 
double  Ra [9] ;/*<*/RMatrix(Ra, a, 2*3. 14159/(200-1)); 

for (int  i=l;i<200;++i) Rotate ( copy (V[i-1] [ j],V[i-l] [ j]+3,V[i] [j])-3,o,Ra);} 
double  R[9] ;/*<-*/RMatrix(R,x, - .5); 

for(int  i=0, j;i<200;++i)for( j=0; j<23;++j)//tilt  top  of  sphere  towards  observer 
Rotate(H[i] [j],o,R),Rotate(V[i] [j],o,R); 

FILE  *f=freopen("sphere.txt", "w",stdout);// . redirect  output  to  a  file 

for(int  i=0;i<200;++i)for(int  j=0; j<23;++j){ 

for (int  k=0;k<3;++k)printf ("%6. 3f , " jH[i] [ j] [k] ); 

for (int  k=0;k<3;++k)printf("%6.3f%s",V[i] [ j] [k], j==23-l&&k==2?"\n" :",");} 
fclose(f ) ; 

}// - YAG  E  NAUT  @GMAI  L .  COM - LAST~UPDATED~02MAY2013 - 


EXAMPLE  2 


#include  <cstdio>// . FILE,freopen(), stdout,printf (),fclose( ) 

ttinclude  <cstdlib>/ / . rand( ) ,  RAND_MAX 

#include  "y_3d_ops .  h"// . RMatrix( ) ,  Rotate ( ) ,  IParameters ( ) ,  Intersect ( ) 

int  main( ){//<================USE  ROTATIONS  AND  INTERSECTIONS  TO  DRAW  A  PINWHEEL 


using  namespace  y3D0ps; 
double  T[4][9]={0,0,0  ,  0,1,0  ,  -.5, .5,0); 
double  T2[4][9]={0,0,0  ,  -.5,. 5,0  ,  -.5,0,0}; 
const  double  z[3]={0,0,l}  ,  o[3]={0,0,0}; 
double  R[9];/*<-*/RMatrix(R,z, 3. 14159/2); 
for(int  i=l;i<4;++i)for(int  j=0; j<3;++j){ 

for (int  k=0;k<3;++k)T2[i] [3*j+k]=T2[i-l] [3*j+k]; /*->*/Rotate(T2[i]+3*j,o, R) ; 
for (int  k=0;k<3;++k)T[i] [3*j+k]=T[i-l] [3*j+k] ; /*->*/Rotate(T[i]+3* j,o, R); } 

FILE  *f=freopen("pinwheel.txt", "w",stdout);// . redirect  output  to  a  file 

for(int  i=0;i<100000;++i){ 

double  S[6]={rand()*2./ RAND_MAX - 1 , r a nd ( ) *2 . / RAND_MAX - 1 , 

1,S[0],S[1],-1>; 
for(int  k=0;k<4;++k){ 

double  t [ 3 ] ;/*<-*/IParameters(t,T[k] ,S); 
double  x[3]; 

if (Intersect (x,t, S)&&t[l]+t [2] <1) 

printf ( "%6. 3f ,%6. 3f ,%6. 3f\n",x[0],x[l] ,x[2] ) ; }} 
fclose(f ) ; 

FILE  *g=f reopen("outline.txt", "w",stdout);// . redirect  output  to  a  file 

for(int  i=0;i<4;++i)for(int  j=0; j<3;++j)for(int  k=0;k<3;++k) 
printf ( "%6.3f%s",T[i] [3*j+k], k==2?"\n" :","); 
for(int  i=0;i<4;++i)for(int  j=0; j<3;++j)for(int  k=0;k<3;++k) 
printf ( "%6 . 3f %s " , T2 [ i ] [ 3* j+k ] , k==2? " \n ":",") ; 
fclose(g) ; 

}// - YAG  E  NAUT  @GMAI  L .  COM - LAST~UPDATED~02MAY2013 - 
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